import Image
import numpy as np
from itertools import product
from math import cos, sin


def concat(p, it):
	if type(it) == list:
		return [p] + it
	return [p, it]


def itereBoucleImbrique(deg, n):
	if deg == 0:
		return []
	if deg == 1:
		return range(n)
	else:
		return [concat(p, it) for it in itereBoucleImbrique(deg-1, n) for p in xrange(n)]


def getEstimationsFromCoef(variables, coefs, degre):
	sizeVariables = variables.shape
	nbLignes = sizeVariables[0]
	nbVariables = sizeVariables[1]
	res = np.zeros((nbLignes, 1))
	idCoef = 0
	for puissances in itereBoucleImbrique(nbVariables, degre+1):
		if sum(puissances) <= degre:
			colUnCoef = np.array([coefs[idCoef, 0] for x in xrange(nbLignes)])
			for i, puissance in enumerate(puissances):
				tmp = (variables[:,i]**puissance).reshape((nbLignes, 1))
				colUnCoef = np.prod((colUnCoef, tmp), axis=0).reshape((nbLignes, 1))
			idCoef += 1
			res = np.sum((res, colUnCoef), axis=0)
	return res


def reglinm(Y, variables, degre):
	sizeVariables = variables.shape
	nbLignes = sizeVariables[0]
	nbVariables = sizeVariables[1]
	Y = np.matrix(Y).reshape((nbLignes,1))
	X = np.ones((nbLignes, 1)) #Voir la description de la regression lineaire multiple de wikipedia pour comprendre a quoi correspond la matrice X
	#On ajoute l'intersection
	first = True
	puissancesBoucleImbrique = itereBoucleImbrique(nbVariables, degre+1)
	del puissancesBoucleImbrique[0]
	for puissances in puissancesBoucleImbrique:
		if sum(puissances) <= degre:
			newCol = np.ones((nbLignes, 1))
			for i, puissance in enumerate(puissances):
				tmp = (variables[:,i]**puissance).reshape((nbLignes, 1))
				newCol = np.prod([newCol, tmp], axis=0)
			X = np.concatenate((X, newCol), axis = 1)
	X = np.matrix(X)
	coefs = X.transpose()*X
	coefs = np.linalg.inv(coefs)
	coefs = coefs * X.transpose()
	coefs = coefs * np.matrix(Y)
	print "Modele :"
	ind = 0
    	for row in product(range(degre+1), repeat=nbVariables):
		if sum(row)<=degre:
			val = ""
			for i, power in enumerate(row):
				val += "coef" + `i` + '**' + `power` + ' * '
			#print val[:-3] + ' : ' + `coefs[ind, 0]`
			ind += 1
	Yestime = getEstimationsFromCoef(variables, coefs, degre)
	SCExplique = np.array(Yestime - Y.mean())
	SCExplique = SCExplique**2
	SCExplique = SCExplique.sum()
	tmp = np.array(Y - Yestime)
	SCResiduelle = np.array(Y - Yestime)
	SCResiduelle = SCResiduelle**2
	SCResiduelle = SCResiduelle.sum()
	SCTotal = np.array(Y - Y.mean())
	SCTotal = SCTotal**2
	SCTotal = SCTotal.sum()
	coefficientDeDetermination = SCExplique/SCTotal
	print "SCExplique :", SCExplique.sum()
	print "SCResiduelle :", SCResiduelle.sum()
	print "SCTotal :", SCTotal.sum()
	print "Coefficient de determination :", coefficientDeDetermination
	return coefs



imageFileName = "texture03nonUnif.jpg"
imageFileName = "texture02lisse.jpg"
image = Image.open(imageFileName)
dataImage = list(image.getdata())
width = image.size[0]
height = image.size[1]

dataReg = []
red = []
green = []
blue = []
for i, (r, g, b) in enumerate(dataImage):
	x = i / width
	y = i % width
	dataReg.append([x, y, cos(x), cos(y), cos(5*x), cos(5*y)]) #, cos(x), cos(y)])
	red.append(r)
	green.append(g)
	blue.append(b)

dataReg = np.array(dataReg)


degre = 3
coefsRed = reglinm(red, dataReg, degre)
coefsGreen = reglinm(green, dataReg, degre)
coefsBlue = reglinm(blue, dataReg, degre)

redEstime = getEstimationsFromCoef(dataReg, coefsRed, degre)
greenEstime = getEstimationsFromCoef(dataReg, coefsGreen, degre)
blueEstime = getEstimationsFromCoef(dataReg, coefsBlue, degre)

dataNewImage = zip(redEstime, greenEstime, blueEstime)
newImage = Image.new(image.mode, image.size) 
newImage.putdata(dataNewImage) 

outImageFileName = imageFileName[:-4] + "Out" + imageFileName[-4:]
newImage.save(outImageFileName)
